Rotational dynamics of CO solvated in small He clusters: a quantum Monte Carlo study 
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The rotational dynamics of CO single molecules solvated in small He clusters (CO@Hejv) has been studied 
using Reptation Quantum Monte Carlo for cluster sizes up to TV = 30. Our results are in good agreement with 
the roto-vibrational features of the infrared spectrum recently determined for this system, and provide a deep 
insight into the relation between the structure of the cluster and its dynamics. Simulations for large TV also 
provide a prediction of the effective moment of inertia of CO in the He nano-droplet regime, which has not been 
measured so far. 
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I. INTRODUCTION 



Thanks to recent progresses in Helium nano-droplet iso- 
lation (HENDI) spectroscopy^, the infrared and microwave 
spectra of small molecules solvated in Hejv clusters are 
now becoming accessible in the small- and intermediate-size 
regimes iSiii^ In the particular case of carbon monoxide, the 
roto-vibrational spectrum of the molecule solvated in small 
He clusters (CO@HeAr) has been recently studied in the size 
range TV = 2 — 20. 5 The infrared spectrum in the 2145 cm -1 
region of the C-O stretch consists of two R(0) transitions 
which smoothly correlate with the a-type (K = 0^0) 
and 6-type (K = 1 «— 0) R(0) lines of the binary complex, 
CO@Hei. The series of o-type transitions — which starts off 
about 7 times stronger for TV = 1 — progressively looses in- 
tensity as TV increases, until it disappears around TV = 7 — 8. 
Around this size, just before it disappears, the o-type line 
seems to split in two. Analogously, around TV = 15 the a-type 
line also seems to split, and the assignment of experimental 
lines becomes uncertain for larger clusters. Elucidating the 
relation existing between the position, number, and intensity 
of the rotational lines and the size and structure of the cluster 
is the principal goal of the present paper. 

Computer simulations of quantum many-body systems 
have also considerably progressed in recent years, allowing 
in some cases to determine the low-lying spectrum of excited 
states. The rotational dynamics of small molecules solvated 
in He clusters and nanodroplets is one of these favorable in- 
stances. The scarcity of low-lying excited states typical of 
superfluid systems makes it possible in this case to extract in- 
formation on the location and intensity of the spectral lines 
from an analysis of the time series generated by quantum 
Monte Carlo random walks^i^ The rotational spectr um o f 
OCS@HeAr has been studied along these lines in Refs. 
Among the many different flavors of quantum Monte Carlo 
available in the literature, we adopt Reptation Quantum Monte 
Cario&ii which we believe presents distinctive advantages in 
the present case and which will be briefly introduced in Sec. 
HP In Sec. milwe present and discuss our results, whereas Sec. 
HVI contains our conclusions. 



II. THEORY, ALGORITHMS, AND TECHNICAL DETAILS 

Virtually all the ground-state quantum simulation meth- 
ods are based on the prior knowledge of some approximate 
wave-function, <£>0i f° r the system under study. In the Vari- 
ational Monte Carlo Method (VMC) one contents oneself 
of this knowledge and the simulation simply aims at cal- 
culating the complicated multi-dimensional integrals which 
are needed to estimate ground-state approximate expecta- 
tion values, (<&o|j4.|<I>o) (here and in the following quantum- 
mechanical operators are indicated with a hat, ~). To this end, 
a random walk in configuration space is generated according 
to the Langevin equation: 



dx = e fn(x) + d£, 



(1) 



where x = {x a } indicates the coordinates of the system, e is 
the step of time discretization, 



f (x) = 2 



G>log($ (x)) 

dx 



(2) 



and d£ is a Gaussian random variable of variance 2e: 
(d£ a d4p) = 2eS a p. Approximate ground-state quantum ex- 
pectation values are then estimated as time averages over the 
random walk, Eq. Q. 

Within Reptation Quantum Monte Carlo (RQMC), exact 
ground-state expectation values and imaginary-time correla- 
tion functions are calculated as appropriate derivatives of the 
pseudo partition function, in the low temperature (large T) 
limit: 



Zo = (*o|e 



-TH 



(3) 



where H is the Hamiltonian of the system. By breaking the 
time T into P intervals of length e = T/P, Eq. Q can be 
given a path-integral representation: 

Z « /$ (x )nf =0 1 (x 4 |e- e5 |x. t+1 )$o(xp)d p + 1 x. (4) 



For the relatively small systems considered here, it is suffi- 
cient to use the primitive approximation to the imaginary-time 
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propagator: 

(s\e- e "\y) oc c ^-y) 2 ne-e[v{, )+ vm/2 + ^ (5) 

where V(x) is the potential energy at point x. The dynamical 
variables of the statistical-mechanical system whose partition 
function is given by Eq. (0} are segments of the VMC random 
walk generated from Eq. ([Q, x(t) of length T, which we call 
reptiles. As the random walk proceeds, the reptile is allowed 
to creep back and forth: new configurations of the reptile are 
accepted or rejected according to a Metropolis test made on 
the integrand of the path-integral representation of Z , Eq. 
0. It can be demonstrated^ that in the large-T limit the 
sample of reptile configurations thus generated is such that 
the sample average of quantities like 



A[x(r)} = I j\(x(r))d, 



(6) 



converges without any systematic bias (but those due to the 
finite values of e and T) to (A) = (* |^|*o), *o being the 
exact ground-state wavefunction of the system. Even more 
interesting is the fact that sample averages of reptile time cor- 
relations, 



Ca(t) = 



1 



T-T 



T-t 



A(x(t'))A(x(t' + T))dT' , (7) 



provide equally unbiased estimates of the corresponding 
quantum correlation functions in imaginary time: 



C A {it) = <tt o |A(it)A(0)|*o) 

= (* |e St 2e- 5t l|* ), 



(8) 



Ca{t) ~ C-Air). The absorption spectrum of a molecule 
solvated in a non polar environment is given by the Fourier 
transform of the autocorrelation function of its electric dipole, 
d: 



cx 2 7 r^|(* |d|*„)| 2 ( 5(i;„-£;o-^) 



e lult {A{t) -d(0))dt, 



(9) 



where ^>o an d are ground- and excited-state wavefunc- 
tions of the system respectively, and Eq and E n the corre- 
sponding energies. The dipole of a linear molecule — such 
as CO — is oriented along its axis, so that the optical ac- 
tivity is essentially determined by the autocorrelation func- 
tion of the molecular orientation versor: c(t) = Cg(i) = 

(^o\e iHi ne~ im n\^o)- We have seen that RQMC gives easy 
access to the analytic continuation to imaginary time of cor- 
relation functions of this kind. From now on, when referring 
to time correlation functions, we will mean reptile time corre- 
lations, i.e. quantum correlation functions in imaginary time. 
Continuation to imaginary time transforms the oscillatory be- 
havior of the real-time correlation function — which is respon- 
sible for the (5-like peaks in its Fourier transform — into a sum 



of decaying exponentials whose decay constants are the exci- 
tation energies, E n — Eq, and whose spectral weights are pro- 
portional to the absorption oscillator strengths, |{\E f o|d|\E' rl )| 2 . 
Dipole selection rules imply that only states with J = 1 can 
be optically excited from the ground state which has J = 0. 
Information on excited states with different angular momenta, 
J, can be easily obtained from the multipole correlation func- 
tions, cj(t), defined as the reptile time correlations of the 
Legendre polynomials: 

cj(r) = (Pj(n(r) -n(0))) 

j \ 



-in 



2J- 



T J2 YJ M {n(r))Y JM (n(0))) (10) 



M=-J 



Both the He-He and the He-CO interactions used here 
are derived from accurate quantum-chemical calculationspiiiA^ 
The CO molecule is allowed to perform translational and ro- 
tational motions, but it is assumed to be rigid. The trial wave- 
function is chosen to be of the Jastrow form: 



$o = exp 



N 



N 



j=i 



(ID 



where r.; is the position of the i-th atom with respect to the cen- 
ter of mass of the molecule, r, = jr.; |, 6i is the angle between 
the molecular axis and r^, and ry is the distance between the 
i-th and the j-th helium atoms. U\ is expressed as a sum of five 
products of radial functions times Legendre polynomials. All 
radial functions (including U2) are optimized independently 
for each cluster size with respect to a total of 27 variational 
parameters. The propagation time is set to T = 1 K^ 1 , with 
a time step of e = 10~ 3 K _1 . The effects of the length of 
the time step and of the projection time have been estimated 
by test simulations performed by halving the former or dou- 
bling the latter. These effects were barely detectable on the 
total energy, and very small on the excitation energies dis- 
cussed below (we estimate that more converged simulations 
would actually improve the already excellent agreement with 
experimentally observed spectra). 

The estimate of excitation energies and spectral weights 
from imaginary-time correlations amounts to performing an 
inverse Laplace transform, a notoriously ill-conditioned op- 
eration which is severely hindered by statistical noised For 
each value of J, we extract the value of the two lowest- 



lying excitation energies, e 



a,b 



-i.e. the two smallest de- 



cay constants in cj (r) — as well as the corresponding spectral 
weights, A J a b , from a fit of c j (r) to a linear combination of 
three decaying exponentials. This fitting procedure does not 
solve in general the problem of obtaining the spectrum from 
a noisy imaginary-time correlation function. However, if we 
know in advance that very few strong peaks, well separated 
in energy, nearly exhaust the entire spectral weight, their po- 
sition and strength can be reliably estimated from this multi- 
exponential fit. In the present study these favorable conditions 
are usually met, although the limitations of the procedure will 
show in some cases, as discussed below. 
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III. RESULTS AND DISCUSSION 

RQMC simulations have been performed for CO@HeA/ 
clusters in the size range N = 1 — 30. In Fig. ^ we report the 
values of the He atomic binding energy, AEn = En-\—En, 
as a function of the cluster size. AEn first increases up to 
TV = 4—5, and it stays roughly constant in the range N = 5—8; 
from this size on AEn starts decreasing, first slowly, then, 
from N = 10 — 11, rapidly down to a minimum at iV = 19. 
For iV > 19 AEn increases again and slowly tends to the 
nanodroplet regime (where it coincides with the bulk chem- 
ical potential, fi = 7.4 K — ) which is however attained for 
much larger cluster sizes than explored hereiii This behav- 
ior can be understood by comparing the shape of the CO-He 
potential energy function, v(r), with the incremental atomic 
density distributions, Apn(t) = Pn(t) — PJV-i(r), where 
Pn is the expectation value of the He density operator: 



N 



(12) 



i=l 



(see Fig. For very small N the atomic binding energy 
is dominated by the He-CO attraction which is strongest in 
a well located atop the oxygen atom. As He atoms fill this 
well, AEn first slightly increases, as a consequence of the at- 
tractive He-He interaction, then, for larger N, the increased 
He-He interaction is counter-balanced by the spill-out of He 
atoms off the main attractive well, until for N x 9 the re- 
duction of the He-CO interaction overcomes the increased at- 
traction and the binding energy starts decreasing steeply. For 
N in the range 10-14 He density accumulates towards the C 
pole, while, around N = 15, the first solvation shell is com- 
pleted and the differential atomic density, ApN, is consider- 
ably more diffuse starting from N — 16. AEn reaches a 
minimum at N = 19. For larger sizes, the trend in the atomic 
binding energy is dominated by the increase of the He-He at- 
traction related to the increase of the cluster size, until it will 
converge to the bulk chemical potential. 
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FIG. 1: Atomic binding energy, AEn = En — En-i, as a function 
of the cluster size in CO @ He at. The horizontal line on the right of 
the figure indicates the chemical potential in bulk 4 He, p = 7.4 K 
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FIG. 2: (color) Upper left panel: He-CO interaction potential. C 
(blue) and O (cyan) atoms are pictured by two circles whose radius 
is the corresponding Van de Waals radius. The other panels picture 
the differential He density, Apjv = Pat — pjv-i for various sizes of 
the CO @ He at cluster. Color convention is rainbow: red to purple in 
order of increasing magnitude. 



In Fig. [5] we report the positions and spectral weights of 
the rotational lines, as functions of the cluster size, N. In the 
size range N = 1 — 9, analysis of the dipole time correlations 
clearly reveals the presence of two peaks, with the weight of 
the higher-energy (6-type) rapidly decreasing by almost a fac- 
tor 2. Note that the sum of the spectral weights of these two 
lines nicely sums to one, indicating that they exhaust all the 
oscillator strength available for optical transitions originating 
from the ground state. For TV between 10 and 12 (shaded area 
in Fig. [5} the situation is less clear. As the weight of the 6-type 
line drops to zero, the statistical noise on its position grows 
enormously. Furthermore the multi-exponential fit introduces 
some ambiguity, as the results are somewhat sensitive to the 
number of terms in the sum. However the important infor- 
mation that one line disappears between N = 10 and 12 is 
clear. For larger N only one relevant line remains, and the ro- 
bustness of the fitting procedure is recovered, with the minor 
exception of the sizes around N = 16, where the minimum 
of the \ 2 appears to be less sharp, possibly correlating with 
the splitting of the line observed in the infrared spectra for 
A^ = 15 (see below). 

In the upper panel of Fig. 0]we compare the rotational struc- 
ture of the observed infrared (vibrational) spectrum- with the 
rotational excitation energies calculated in this work. Exper- 
imental data are referred to the center, vq, of the vibrational 
band for N = (CO monomer). In order to better compare 
our predictions with experiments, we have corrected the for- 
mer with an estimate of the vibrational shift, Afo, i.e. the dis- 
placement of the vibrational band origin as a function of the 
number of He atoms. The vibrational shift can be calculated as 
the difference in the total energy of the cluster obtained with 
two slightly different potentials^ vqo and vxi, representing 
the interaction of a He atom with the CO molecule in its vibra- 
tional ground state and first excited states, respectively. Our 
estimate of the vibrational shift as a function of the cluster size 
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FIG. 3: Upper panel: position of the rotational lines of CO@Hejv 
as obtained from RQMC simulations as a function of the cluster size, 
N. a-type lines are indicated with triangles, o-type lines with dia- 
monds. Lower panel: spectral weights of the lines reported in the 
lower panel; the continuous line near the upper border of the figure 
corresponds to the sum of the spectral weights. 



is reported in the lower panel of Fig. [4] Since the evaluation 
of a small difference between two large energies is computa- 
tionally demanding for large clusters, A^o has been evaluated 
perturbatively with respect to the difference t>oo — i>ii We 
have used the vibrational shift calculated in Ref. [3 after ver- 
ifying on small clusters that the perturbative treatment is reli- 
able. The agreement between our results and experiments is 
remarkable. Some of the features of the observed spectrum, 
however, call for a deeper understanding and theoretical in- 
vestigation. Two questions, in particular, naturally arise. Why 
two peaks are observed in the small-size regime, and what de- 
termines the disappearance of one of them at N = 8? What 
determines the split of the higher-frequency (6-type) line at 
N = 7 and of the lower-frequency (a-type) one at N = 15? 

The existence of two lines for small N is likely due to 
a larger asymmetry of the cluster in this regime. If the 
CO@HeAr complex is described as a rigid rotor, in fact, one 
would have one rotational line originating from a J = 
ground state if the complex has cylindrical symmetry, while 
this line would double if some of the atomic density accumu- 
lates in a longitudinal protrusion. The inertia of the complex 
would in this case be larger for a rotation about an axis perpen- 
dicular to a plane containing the protrusion (end-over-end ro- 
tation) than about an axis lying on such a plane. Given that the 
He density in the ground state of CO@HeAr is cylindrically 
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FIG. 4: Upper panel: positions of the infrared lines of CO @ He n 
as observed experimentally (Ref. empty circles) and as estimated 
from the present simulations and corrected by the estimated vibra- 
tional shift (solid diamonds, see text). Lower panel: vibrational shift 
of the lines, as estimated in the present work (triangles) and in Ref. 
FLU (diamonds and dashed line). 



symmetric, any departure from this symmetry can only show 
up in higher correlation functions. The situation is conceptu- 
ally similar to that of a fluid whose density is homogeneous 
and whose structure at the atomic scale is reflected in the pair 
correlation function. Analogously, we define an atomic an- 
gular correlation function, C(<p), as the probability of finding 
two He atoms which form a dihedral angle <fr with respect to 
the molecular axis: 



C{<f>) 



(N- 1) 



(13) 



In the upper panel of Fig. [5]we show C (4>), for different clus- 
ter sizes. The depletion of C for larger than 7r/2, clearly 
visible for N = 3 (green circles), indicates a tendency of the 
He atoms to cluster on a same side of the molecular axis. For 
larger clusters, however, this effect weakens to the extent that 
it becomes difficult to disentangle from the structural infor- 
mation related to the He-He interaction (the dimple at small 
and the subsequent maximum around <fi = 0.3 — 0.4 if). A 
more sensitive measure of the propensity of He atoms to clus- 
ter on a side of the molecule is given by the integral of C(</>) 
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FIG. 5: (color) Upper panel: probability density of finding two He 
atoms which form a dihedral angle <j> with respect to the molecular 
axis; the probability is normalized to 1. Results pertain to clusters 
with TV = 3 (green circles), N = 6 (blue triangles) and TV = 13 (red 
diamonds). Lower panel: the integrated probability density, defined 
in Eq. J 141 . as a function of the cluster size. 

from to §, 

M= f* C(4>)d4>-\. (14) 

In the lower panel of Fig. [5] we display M as a function of 
the cluster size, TV: one sees that M decreases with TV and 
reaches a minimum at TV ~ 14. This is the size at which 
the first solvation shell is completed, and the cluster asym- 
metry increases again when the second shell starts to build. 
The rotational spectrum of the solvated molecule, however, 
is insensitive to this asymmetry for clusters of this and larger 
sizes because the motion of He atoms in the second and outer 
solvation shells is decoupled from that in the first and from 
molecular rotation. The existence of a longitudinal asymme- 
try is a necessary condition for the doubling of the rotational 
line. Whether or not this condition is also sufficient depends 
on the dynamics: if quantum fluctuations make the motion 
of the protrusion around the molecular axis fast with respect 
to the molecular rotation, then the asymmetry is effectively 
washed out. The existence of two lines in the rotational spec- 
trum of the molecule implies therefore that an asymmetry in 
the classical distribution of He atoms around the molecular 
axis exists; that the molecular inertia is sensitive to this asym- 
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FIG. 6: Diamonds: CO rotational frequencies in CO@Hejv as func- 
tions of the cluster size, N (same as in Fig. |3j- Dots: frequency 
of the lowest mode appearing in the spectral analysis of the angular 
He-He correlation function (see Eq. i 151 V 

metry (the protrusion can be 'dragged' along the molecular 
rotation); and that the motion of this protrusion around the 
molecular axis is not adiabatically decoupled from the molec- 
ular rotation. 

In order to better characterize the motion of He atoms 
around the molecule and the coupling of this motion to molec- 
ular rotation, we examine the imaginary-time correlations of 
the versor, u, of the He center of mass, Tcm> relative to the 
molecular center of mass: 

C u (r) = (u(r) • u(0)>. (15) 

For the binary complex, He-CO, vqm coincides with the po- 
sition of the helium atom, and we expect its angular dynamics 
to be strongly correlated to the molecular rotation, at least in 
the end-over-end mode. In Fig. [6]we report the frequency of 
the slowest mode appearing in the spectral analysis of C u (t), 
e u , as a function of TV, and compare it with the corresponding 
frequencies of the molecular rotation. We see that for cluster 
sizes up to TV = 9 — 10, e u is degenerate with the a-type fre- 
quency in the molecular rotational spectrum, with a spectral 
weight which passes from A u w 1 for TV = 1 to A u « 0.7 for 
TV = 10. These findings are a manifestation of the fact that He 
atoms are dragged along the slowest, end-over-end, rotation 
of the solvated molecule, and that the effect of this dragging 
decreases when more He states with J = 1 become available 
and subtract spectral weight to the slowest mode. For TV > 10, 
e u further increases and departs from e a , indicating an effec- 
tive decoupling of the two kinds of motion. In this regime, 
the effective rotational constant B of the solvated molecule is 
almost independent of the cluster size. Free molecular rota- 
tion with an increased moment of inertia with respect to the 
gas phase is the typical signature of superfluid behavior in He 
nanodroplets. Extrapolating the result obtained for TV up to 
30 to the nanodroplet limit, we predict a renormalization fac- 
tor of the B value of 0.78. The lowest atomic mode, e u , slows 
down again for TV = 15. This is due, however to the slow He 
motion in the second solvation shell which hardly affects the 
rotation of the solvated molecule. Although the resolution that 



6 



can be achieved with our simulations is not sufficient to detect 
the doubling of the a and 6 lines which is experimentally ob- 
served for TV = 15 and TV = 7 respectively, it is interesting 
to notice that the former occurs in correspondence with the 
crossing between e u and e a , possibly due to the resonant in- 
teraction between the two modes. It is tempting to assume 
that a similar mechanism may be responsible for the doubling 
of the b line at TV = 7, involving however higher-energy He 
states. A deeper study of the He dynamics would clarify this 
point. 

IV. CONCLUSIONS 

Computer simulations of quantum many-body systems 
have reached such a degree of sophistication and reliability 
that in some cases they can be used to provide information, 
complementary to that which can be obtained in the labora- 
tory, on the dynamical processes probed spectroscopically. 

In the case of small polar molecules solvated in He clusters, 
for instance, the calculation of the time autocorrelation of the 
molecular dipole (which is the quantity directly coupled to the 
experimental probe) allows to reproduce rather accurately the 
roto-vibrational excitation energies which are now becoming 
experimentally accessible for small clusters (TV = 1 — 20). 
Even more importantly, computer simulations give direct ac- 
cess to quantitities and features (such as, e.g., static and dy- 
namic properties of the He matrix) which are not accessible to 
the experiment, and whose knowledge provides the basis for 
understanding the relation between structure and dynamics in 
these confined boson systems. 

In the specific case of CO@HeAr, which is the subject 
of the present study and of a recent infrared spectroscopy 
experiment^ 5 the presence of two spectral lines — a-type and 
6-type, evolving respectively from the end-over-end and from 
the free-molecule rotations of the binary complex — is related 



to the propensity of the He atoms to cluster on a same side 
of the molecular axis, which we measure by an angular pair 
distribution function: as more He atoms progressively fill the 
first solvation shell, their clustering propensity weakens; the 
CO impurity gets more isotropically coated, looses a preferred 
axis for the free-molecule mode, and the &-type line disap- 
pears. 

The time autocorrelation of the versor of the He center of 
mass provides dynamical information on the He atoms in ex- 
cited states with J = 1. We find a substantial spectral weight 
on a He mode whose energy, e u , is degenerate with the a-type 
line for TV up to about ten. This indicates that some of the He 
density is dragged along by the molecular rotation — in other 
words, part of the angular momentum in the cluster mode in- 
volving molecular rotation is carried by the He atoms. We also 
find that for larger clusters the molecular rotation effectively 
decouples from this He mode, and its energy e a becomes es- 
sentially independent of the number of He atoms. Based on 
the nearly constant value of e a in the range of TV between 15 
and 30, well beyond completion of the first solvation shell, 
we predict the effective rotational constant in the nanodroplet 
limit to be smaller by a factor 0.78 than its gas phase value. 
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